#define _EPS 0.01
#define NB_K1  2
#define P_K 1
//#define MAXITK 80
#include "../Particle.hpp"
#include "Kernel_Probit.hpp"
template<class Resample>
class Particle_Probit{

public:
	Particle_Probit(Kernel_Probit<ProbitTemp> *K,Resample *R,ProbitTemp *D,int M, Distribution::Distribution *F){
		_p=(*D).Get_p();
		mat Xt(M,_p);//p le nombre de parametre a obtenir de Density
		W=new double[M];
		X=Xt;
		_K=K;
		_R=R;
		_F=F;
		_Z=0;
		_D=D;
		kk=-1;
		_thres=_C*M;
		_M=M;
		Normalize();
		y=X.t();
		_n=(*D).Get_n();
		_K->Set_n(M);
	//	cout << "///" << _n << "///";
	}
	Particle_Probit(){
		cout << "Hidee2!";
	}
	Particle_Probit(const Particle_Probit& X)
	{
		cout << "Hidee!";
	}
	~Particle_Probit(){
	//	delete[] W;
	}
	void Filter()
	{
		this->Init();
		int i=0;
		double p=0; 
		//for(int i=1;i<(_n);i++)
	//	i{
		double b=(*_D).Get_bn();
		_temp=0;
		Phiv.push_back(_temp);
		b+=-_temp;
		//b=Cooling(b,i);
		//_temp=StepLength(b);
		//Phiv.push_back(_temp);
		//b+=_temp;
		cout << "b: " << b << " p";
		_D->Set_xi(b);
		y=X.t();
		_D->Weight(y,W,i);//le p designera comment on avance 
		double ess=Ess_W(W,1);
		cout << " ESS: " << ess;
		X=y.t();
		p=_D->Get_bn();
		while((p-1)!=0){
			i++;
		//	cout << i;
			cout << "step: " << i << "\n";
			this->Step(i);
			p=_D->Get_bn();
			cout << "the p "<< p-1;
			if(p==1){cout << "Hello";}
		}
		
	}
	void Init()
	{
		//pior
		X=_F->scrambled(_M);	
	}
	void Step(int i)
	{
		//cout << "Ess:" << Ess << "\n";
			double Ess=Ess_W(W,0);
	//cout << "\\";
		if(Ess< _thres)
		{
			mat v=growingvect(_M);
			double *w=new double[_M];
			for(int j=0;j<_M;j++){
				 w[j]=exp(W[j]);
			//	cout << w[j] << "\n";
			}
			(*_R)(&_M,w,v);
			delete[] w;
			//cout << v << "\n";
			Arangemat(X,v);
			Normalize();//poids fixé a 1
				
			
			
			for(int o=0;o<NB_K1;o++){	
				cout << "avant " << mean(X,0); 
				_K->Move(&(X),i);
				cout << "après " << mean(X,0); 
				
				cout << i;
			}
		}
				
			y=X.t();
			//cout << X << "////";	
			double b=(*_D).Get_bn();
			_temp=StepLength(b);
			Phiv.push_back(_temp);
			b+=-_temp;

			if(b<1){b=1;}
			cout << "b: " << b << "p";
			_D->Set_xi(b);

		//cout<< "before:" << W[0] ;
		//cout<< "after:" << W[0] ;
		//dans le cas de tempering peut etre determiner
		_D->Weight(y,W,i);//le p designera comment on avance 
	//	cout << "after:" << y.col(1);
		X=y.t();
		double ess=Ess_W(W,1);
		cout << " ESS: " << ess;
		//cout << " ESS: " << Ess_W(W,1);
	}

	double Ess_W(double *w, int tf)//renormalisation???
	{
		double *foo= new double[_M];
		//double *foo2= new double[_M];
		double sum=w[0];
		double sum2=exp(w[0]);
		for(int i=1;i<_M;i++)
		{
			double t=log_add(w[i],sum);
			sum=t;
			//sum2+=exp(w[i]);
		}
		//cout << "sum " << sum-log(sum2);

		if(tf){
			_Z+=sum-log(_M);
		}
		//cout << "the sum "<< sum;
		
		for(int j=0;j<_M;j++)
		{
			foo[j]=(w[j]-sum);
		//	foo2[j]=exp(w[j])/sum2;
			//cout << exp(foo[j])-exp(w[j])/sum2 << "\n";
		}

		sum=2*foo[0];
		sum2=foo[0]*foo[0];
		for(int i=1;i<_M;i++)
		{
			double t=log_add(sum,2*foo[i]);
			sum=t;
		//	sum2+=foo2[i]*foo2[i];
		}
		delete[] foo;
		//delete[] foo2;
		//cout << "logESS "<< -sum;
		return exp(-sum);		
	}

	double Cooling(double b, int i)
	{
		double s=0;
		if(i<10)
		{
			s=0.0005;
		}else if(i<40){
			s=0.001;
		}else{
			s=0.005;
		}
		return (double)1/((double)1/b+s);
	}
	double StepLength(double phi)
	{
		double u=phi-1;
		double l=0;
		double thres=0.000001;
		double e=2;
		double alpha=0.001;
		double eta=0;	
		int Bool=1;
		while((e>thres)&(l<phi-1))
		{
			mat y=X.t();
			//cout << "eta"; 
			//for(int i=0;i<_M;i++)
			//{
			 //      	cout << [i] << "\n";
			//}
			eta=Ess_W(_D->Eta_star(alpha,W,y,Bool),0)/_M;
			if(eta>_U)
			{
				l=alpha;
				alpha=(double)(alpha+u)/2;
			}else{
				u=alpha;
				alpha=(double)(alpha+l)/2;

			}
			Bool=0;	
			e=abs<double>(u-l);
			cout << "\\\\\\ "<< eta <<" \\\\\\\\\n"; 
		}
		if(1>phi-alpha)
		{
			alpha=phi-1;
		}

		return alpha;
	}
		
	mat Get_theta(void){return X;}
	void Normalize(void){
		for(int i=0;i<_M;i++){
			W[i]=0;
		}
	}
	mat   Get_W(void){
		double sum=0;
		mat C(_M,1);
		for(int i=0;i<_M;i++)
		{
			sum+=W[i];
		}	
		for(int i=0;i<_M;i++)
		{
			C(i,0)=(W[i]/sum)*_M;
		}
//		cout << C;
		return C;
	}
	double Get_Z(void){
		cout << "ev: " << _Z;       
		return _Z;}
	mat Get_Phiv(void)
	{
		int n=Phiv.size();
		mat res(n,1);
		for(int i=0;i<n;i++)
		{
			res(i,0)=Phiv[i];
		}
		return res; 
	}
	mat Get_Acc(void)
	{
		int n=Acc.size();
		mat res(n,1);
		for(int i=0;i<n;i++)
		{
			res(i,0)=Acc[i];
		}
		return res; 
	}

private:
	mat X;//value particle
	double *W;//weigths
	double Ess;
	int _M;
	double _thres;
	double _Z;
	Kernel_Probit<ProbitTemp> *_K;
	Resample *_R;
	ProbitTemp *_D;	
	Distribution::Distribution *_F;
	mat y;
	double _temp;
	int _n;
	int _p;
	int kk;
	vector<double> Phiv;
	vector<double> Acc;
};
